function FVSWE1Dv2 
% File name: FVSWE1Dv2.m 
% This m-file is associated with the free text book: 'Introductory Finite  
% Volume Methods for PDEs' which may be downloaded from bookboon.com 
% when published. 
% 
% Notes:  
% 1. This prog is version 2 of a sequence of programs which gradually  
% build up to coding a FV solver on a general 2D structured mesh.   
%  
% 2. This code is based on FVSWE1Dv1.m by Clive Mingham.  The changes are 
% that mesh cells are given by x and y coordinates of their vertices 
% and cell areas and side vectors are stored in arrays and calculated 
% for general quadrilateral cells although only left and right side 
% vectors are found in this 1D case. 
% 
% Description:  Solves the system of 1D PDEs, dU/dt + div(H) = 0, 
% using a finite volume scheme on a single row of uniform Cartesian cells  
% in the x direction.   
% 
% This equation corresponds to the finite volume form: 
% 
% doubleintegral(dU/dt dA) + lineintegral(H.n ds) = 0. 
% 
% In this case the equations are the 1D shallow water equations (SWE) so:  
% Water depth: h=h(t,x) [m]. 
% Water velocity: v = vx i + 0 j where vx=vx(t,x) [m/s]. 
% Acceleration due to gravity: g=9.81 [m/s^2]. 
% Geopotential: phi = gh.  
% Matrix of conserved variables: U = [ phi 
%                                     phi vx]. 
% Flux density: H = [       phi v           i + [0   j 
%                    phi vx^2 + 0.5*phi^2 ]      0] 
% 
% subfunctions: freadmesh, fcellareas, fsidevectors, finitialU,  
%               finterfaceflux, finsertghostcells, fcalcdt,  
%               fsolver, fremoveghostcells, fplotresults, fgetHifromUi. 
% 
% Verification:  results look exactly the same as for FVSWE1Dv1.m 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% 
clc; clf; clear all; % clear everything 
runtime=1.2;  % runtime [s] 
t=0;          % current time [s] 
timelevel=0;  % current time level 
[x,y,xcen,ycen]=freadmesh; % gets coordinates of the computational mesh 
                           % and cell centres 
A=fcellareas(x,y);   % compute and store areas for each cell.  Includes 
                     % left and right ghost cells. 
S=fsidevectors(x,y); % compute and store cell side vectors in an 
                     % (NI+2)x2x2 array which contains left and  
                     %  right ghost cells. 
% 
U=finitialU(xcen); % Put initial cell centre values of U in a 2xNI array 
% 
Uinitial=U;  % Store initial profile for plotting if required 
U=finsertghostcells(U); % Append ghost cells to each end of the domain.  
                        % This means that indices of computational cell  
                        % start at i=2 and end at i = NI+1. 
Unext=zeros(size(U)); % Creates an array for the updated U values  
[d,N]=size(U);        % d is the number of PDEs (=rows of U) 
NI=N-2; % number of computational cells in the x direction 
%% 
disp('time marching starts')  
while(t<runtime)    
 timelevel=timelevel+1  
 U=fbcs(U);   % Insert boundary conditions into ghost cells at each end 
              % of the domain. U becomes 2x(NI+2) so each computational  
              % cell is indexed i=2 to NI+1. 
 dt=fcalcdt(A,S,U); % Finds stable time step before each iteration. 
 t=t+dt;  % updates time   
 if (t>runtime) % ensures program stops at exactly t=runtime [s] 
     dt=dt-(t-runtime); 
     t=runtime; 
 end 
 t   % Outputs current time to the Command Window 
 dt 
 Unext=fsolver(U,A,dt,S); % implements a time step of the solver 
% 
 U(:,2:NI+1)=Unext(:,2:NI+1);  % update computational U values prior to 
                               % the next time iteration 
end  % of while loop 
disp('1D finite volume scheme ended at time t') 
t 
U=fremoveghostcells(U); 
fplotresults(xcen,U) % plot results 
disp('finished main program, see figures') 
end % FVSWE1Dv2 
%%%%%%%%%%%%%%%%%%%%%%%%%%% subfunctions  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [x,y,xcen,ycen]=freadmesh 
% Verification:  verified. 
% The mesh has NI uniform quadrilateral cells in the x direction and  
% NJ=1 cell in the j direction (for a 1D problem). 
% The x and y coordinates of the lower left hand corner of cell (i,j) are 
% held in arrays x and y respectively which are both (NI+1)x(NJ+1).  In 
% this way the 4 vertices of cell (i,j) are (x(i),y(j)), (x(i+1),y(j)), 
% (x(i+1),y(j+1)) and (x(i),y(j+1)).   
% The mesh starts at the origin so x(1,1)=0=y(1,1). 
%  
NI=100;  % number of cells in the i direction 
NJ=1;    % number of cells in the j direction (=1 for a 1D problem) 
x=zeros(NI+1,NJ+1); % allocate correct sized array for cell vertices 
y=zeros(NI+1,NJ+1); % allocate correct sized array for cell vertices 
xcen=zeros(NI,NJ); % allocate correct sized array for cell centres 
ycen=zeros(NI,NJ); % allocate correct sized array for cell centres 
% 
dx=1.0; % cell length in x direction 
dy=2.0; % cell length in y direction (arbitrary) 
% 
x(1,1)=0; % first row starting point for x 
y(1,1)=0; % first row starting point for x 
x(1,2)=0; % second row starting point for x 
y(1,2)=dy;  % second row starting point for y 
% create uniform mesh 
for i=2:NI+1 
      x(i,1)=x(1,1)+(i-1)*dx; 
      y(i,1)=y(1,1); 
      x(i,2)=x(1,2)+(i-1)*dx; 
      y(i,2)=y(1,2); 
end % of i loop 
% find cell centres by averaging coordinates of vertices  
% (this works for a general structured mesh) 
for i=1:NI 
     xcen(i,1)=(x(i,1)+x(i+1,1)+x(i+1,2)+x(i,2))/4; 
     ycen(i,1)=(y(i,1)+y(i+1,1)+y(i+1,2)+y(i,2))/4; 
end % of i loop 
end % freadmesh 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function A=fcellareas(x,y) 
% General way to find quadrilateral cell areas. 
% Defines A to include left and right ghost cells. 
 [m,n]=size(x); 
 NI=m-1; % number of computational cells in the x direction 
 NJ=n-1; % number of computational cells in the y direction 
 A=zeros(NI+2);  % allocate correct size array including ghost cells 
 % inclusion of a left ghost cell means that indices of computational cells 
 % go from i=2 to i=NI+1 
 for i=1:NI 
     for j=1:NJ % j goes from 1 to 1 
         % write the diagonal vectors as 3D vectorrs for the cross product 
         d1=[x(i+1,j+1)-x(i,j),y(i+1,j+1)-y(i,j),0];  
         d2=[x(i,j+1)-x(i+1,j),y(i,j+1)-y(i+1,j),0]; 
         d=cross(d1,d2);    
         A(i+1)=0.5*norm(d); % note shift of i index by 1 in A 
     end % of j loop 
 end % of i loop 
end % of fcellareas 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
function S=fsidevectors(x,y) 
% Verification:  
% 
% Finds the left and right side vectors for each cell and stores in array 
% S.  The first component of S holds the cell index.  The second component 
% of S indicates the right (entry = 1) or left (entry = 2) side vector. 
% The third component of S indicates the i (entry = 1) or j (entry = 2)  
% component of the side vector. 
% 
% If a i + b j is the vector representing the cell side then the  
% corresponding side vector is b i - a j. 
% Notes: 
% 1. S contains left and right ghost cells. 
% 
% 2. vectors are 1 by 2 arrays so care must be taken because side 
%    vectors are stored in higher dimensional arrays. 
% 
[m,n]=size(x); 
NI=m-1; % number of computational cells in the x direction 
NJ=n-1; % number of computational cells in the y direction 
S=zeros(NI+2,2,2); % includes ghost cells 
% ghost cells shift the i index by 1 so computational indices 
% go from 2 to NI+1 
for i=1:NI 
    for j=1:NJ        
        rightcellside=[x(i+1,j+1)-x(i+1,j),y(i+1,j+1)-y(i+1,j)]; 
        S(i+1,1,1)=rightcellside(2); % note shift of i index 
        S(i+1,1,2)=-rightcellside(1);        
        % 
        leftcellside=[x(i,j)-x(i,j+1),y(i,j)-y(i,j+1)]; 
        S(i+1,2,1)=leftcellside(2); 
        S(i+1,2,2)=-leftcellside(1);     
    end % of j loop 
end % of i loop 
end % fsidevectors 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function U=finitialU(xcen) 
% This sets up a dam break problem. 
% Verification:  correct. 
% Inserts initial U values in computational cells (U is 2xNI) 
NI=length(xcen); % number of computational cells in x direction 
% 
d=2;  % number of PDEs 
U=zeros(d,NI); % create correct sized array (without ghost cells) 
dx=xcen(2)-xcen(1); 
lengthofdomain=NI*dx; 
xmidpoint=lengthofdomain/2; 
g=9.81; 
% 
% 1D Dam break.  Left and right states are either side  
% of the centre of computational domain. 
% Left states 
hleft=10.0; % left water depth [m] 
vleft=[0,0]; % left water velocity vector [m/s] 
phileft=hleft*g; % left geopotential 
% Right states 
hright=5.0; % right water depth [m] 
vright=[0,0]; % right water velocity vector [m/s] 
phiright=hright*g; % right geopotential 
% 
% Fill in U values for each computational cell 
for i=1:NI 
    if (xcen(i)<=xmidpoint) 
        U(1,i)=phileft; 
        U(2,i)=phileft*vleft(1); 
    else 
        U(1,i)=phiright; 
        U(2,i)=phiright*vright(1); 
    end % of if 
end % of i loop 
% 
end % finitialU 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
function U=finsertghostcells(U) 
% Verification:  verified. 
% Inserts ghost cells at left and right ends of the computational domain.  
% Hence computational indices go from i=2 to NI+1. 
% Assumes that U has 2 rows (i.e. 2 PDEs) which is correct for the 1D SWE. 
%  
 [d,NI]=size(U); % d is the number of equations (= rows of U). 
 leftghost=[0 
            0]; 
% 
 rightghost=[0 
             0]; 
% 
 U=[leftghost U rightghost]; 
end % finsertghostcells 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
function U=fbcs(U) 
% Verification:  verified. 
% Inserts boundary conditions in to the ghost cells at left and right ends  
% of the domain.  Note that ghost cells have already been appended to U and 
% U is initially zero in each ghost cell. 
%  
% The following applies zero gradient boundary conditions everywhere. 
 [d,N]=size(U); % d is the number of equations (= rows of U). 
% Left 
 U(:,1)=U(:,2); % copies the neighbouring interior values to the ghost cells 
% Right 
 U(:,N)=U(:,N-1); % copies the neighbouring interior values to the ghost cells 
end % fbcs 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function Unext=fsolver(U,A,dt,S) 
% Implements the finite volume scheme for 1 time step 
[d,N]=size(U); % includes ghost cells 
NI=N-2; % number of computational cells 
for i=2:NI+1 % loop over computational cells 
     SRight=[S(i,1,1),S(i,1,2)]; % right side vector 
     SLeft=[S(i,2,1),S(i,2,2)]; 
     [IHRight,IHLeft]=finterfaceflux(U,i); % finds the left and right  
     %                                       interface fluxes for cell i 
     % Update U to Unext 
     for k=1:d 
     % compute total flux out of cell i for equation k 
       IHR=[IHRight(k,1),IHRight(k,2)]; % right flux vector 
       IHL=[IHLeft(k,1),IHLeft(k,2)];   % left flux vector 
       totalfluxout=dot(IHR,SRight)+dot(IHL,SLeft); 
       % Uses Lax-Friedrichs 
       Unext(k,i)=0.5*(U(k,i-1)+U(k,i+1))-(dt/A(i))*totalfluxout;  
     end % k loop 
 end % of i loop 
 end % fsolver 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [IHRight,IHLeft]=finterfaceflux(U,i) 
% Interface H values are based on linear interpolation of neighbouring cell  
% centre data which is the same as averaging on a uniform mesh.   
% These may be found either by finding the average value of U at the cell  
% at the interface and then constructing H or by direct averaging of H  
% components. Since the PDEs are non-linear these approaches are not the same. 
% The latter appoach is analogous to the Lax-Friedrichs finite  
% difference method and is the better method as the former method produces 
% unwanted oscillations. 
% 
% Interface fluxes are denoted IHRight, IHLeft. 
% For the 1D shallow water equations H is: 
% [     phi*vx          i  +  [ 0     j 
%    phi*vx^2+phi^2/2 ]         0 ]     
  
[d,m]=size(U);   % d=number of rows of U (= number of PDEs = 2) 
NI=m-2;  % number of computational cells 
IHRight=zeros(d,2); % array to store the right interface flux for a cell.  
                    % IHRight(,1) is the i component of H and 
                    % IHRight(,2) is the j component of H which is zero. 
IHLeft=zeros(d,2); 
% 
% Find interface values of H directly by averaging neighbouring components. 
% Cell i 
Hi=fgetHifromUi(U,i); 
% Cell i-1 
Himinus1=fgetHifromUi(U,i-1); 
% Cell i plus 1 
Hiplus1=fgetHifromUi(U,i+1); 
% 
% Right interface is the average of H components in cells i and i+1. 
IHRight=(Hiplus1+Hi)/2; 
% Left interface is the average of H components in cells i-1 and i. 
IHLeft=(Himinus1+Hi)/2; 
end % finterfaceflux 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
function outarray=fremoveghostcells(inarray) 
% Verification:  Verified 
% Removes left and right ghost cells. 
 [d,m]=size(inarray); 
 outarray=inarray(:,2:m-1); 
end % fremoveghostcells 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function dt=fcalcdt(A,S,U) 
% Verification:  not verified. 
% Finds allowable time step dt using heuristic formula. 
% 0<F<1 is a tunable safety factor. 
 [nrows,M]=size(U);  % U now enlarged to include ghost cells  
 dt=9999999;  % large value of dt to start 
 F=0.9; % safety factor 
 for i=2:M-1 % indices go over computational cells  
     phi=U(1,i); 
     cel=sqrt(phi); % wave celerity 
     vx=U(2,i)/phi; % flow speed in x direction 
     v=[vx,0]; % flow velocity 
     SRight=[S(i,1,1),S(i,1,2)]; 
     a=dot(v,SRight); 
     b=cel*norm(SRight); 
     % 
     speed=max(abs(a+b),abs(a-b));  % max abs wave speed 
     localdt=A(i)/speed; 
     % take smallest value of dt 
     if (localdt<dt) 
        dt=localdt; 
     end 
  end % i loop 
 dt=F*dt; 
end  % fcalcdt 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function fplotresults(xcen,U) 
% Plots results on a structured mesh - even works for non-Cartesian! 
% 
  g=9.81;  % acceleration due to gravity 
  waterdepth=U(1,:)/g; 
  waterspeed=U(2,:)./U(1,:); 
  subplot(2,1,1), plot(xcen,waterdepth) 
  title('FVSWE1Dv2.m: numerical solution for water depth') 
  xlabel('x [m]') 
  ylabel('water depth [m]') 
  % 
  subplot(2,1,2), plot(xcen,waterspeed) 
  title('FVSWE1Dv2.m: numerical solution for x component of water velocity') 
  xlabel('x [m]') 
  ylabel('speed [m/s]') 
% 
end % fplotresults 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function H=fgetHifromUi(U,i) 
% Constructs flux vector H at centre of cell i from U(i)  
% U = [  phi 
%       fi*vx ] 
% 
% H=[     phi*vx    ,  0          Each row of H is a 2D vector  
%    phi*vx^2+fi^2/2,  0 ] 
% 
[d,n]=size(U); % d is the number of PDEs (= rows of U) 
H=zeros(d,2);  % flux density vectors 
% 
phi=U(1,i); % geopotential 
vx=U(2,i)/phi; % velocity component in x direction 
% 
H(1,1)=U(2,i); % i component of first row of H 
H(1,2)=0;      % j component of first row of H 
% 
H(2,1)=phi*vx*vx+0.5*phi*phi; % i component of second row of H 
H(2,2)=0;                     % j component of second row of H 
% 
end % fgetHifromUi 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 

